#function for generating test data table of pairwise distances
#minimizes computation by slimming calculations via a grid
# slower than pdist for small sets (up to 2000 ), otherwise wins hands down
pdist.fast<-function(X, Y, maximum, space.size){
	if (min(X)<0) X<-X-min(X) # move negative coordinates to positive land
	if (min(Y)<0) Y<-Y-min(Y)
	lth<-space.size/maximum
	if (lth%%1!=0) {
		print("Error: maximum must divide into space.size perfectly")
		return(NULL)	
	}
	out<-vector("list", length=length(X)) #list to take neightbour.ID and distance
	points<-vector("list", length=lth^2)#matrix of lists of point IDs
	gX<-X%/%maximum+1 #collapse to grid refs
	gY<-Y%/%maximum+1
	neigh<-neighbours.init(space.size, maximum)
	for (i in 1:length(X)){ #throw point IDs into grid cell list
		temp<-flatmat(gX[i], gY[i], lth)
		points[[temp]]<-c(points[[temp]], i)
	}
	for (i in 1:length(X)){
		nb<-neigh[,,gX[i], gY[i]] #find neighbouring cells
		nb<-subset(nb, is.na(apply(nb,1,sum))==F)
		temp<-flatmat(nb[,1], nb[,2], lth) #find relevant points
		snk.ID<-unlist(points[temp])
		dists<-sqrt((X[i]-X[snk.ID])^2+(Y[i]-Y[snk.ID])^2)
		temp<-cbind(snk.ID, dists)
		temp<-subset(temp, temp[,"dists"]<=maximum)
		out[[i]]<-temp
	}	
	out
}

#takes x, y, indices of a matrix and places them into a vector
flatmat<-function(x, y, size){ #size is size of one side of matrix
	(y-1)*size+x
}


neighbours.init<-function(space.size, cell.size){
	lth<-space.size/cell.size
	neigh<-array(0, dim=c(9,2,lth,lth))	
	for (x in 1:lth){
		for (y in 1:lth){
			out<-matrix(nrow=9, ncol=2)
			ifelse((x-1)==0, X<-NA, X<-x-1) #correct for indices going to zero
			ifelse((y-1)==0, Y<-NA, Y<-y-1)
			ifelse((y+1)>lth, Yp<-NA, Yp<-y+1) # correct for indices going to > space
			ifelse((x+1)>lth, Xp<-NA, Xp<-x+1)
			# Assign neighbours
			out[1,]<-c(x, y) #record target cell
			out[2,]<-c(Xp, y)
			out[3,]<-c(X, y)
			out[4,]<-c(x, Y)
			out[5,]<-c(Xp, Y)	
			out[6,]<-c(X, Y)
			out[7,]<-c(x, Yp)
			out[8,]<-c(Xp, Yp)
			out[9,]<-c(X, Yp)
			neigh[,,x,y]<-out
		}	
	}
	neigh	#return the array
}

setwd("/home1/30/jc227089/Toad_point_process_model/DONE_ART")
d<-read.csv("pop_tabl.csv")
pairs.art<-pdist.fast(d$X, d$Y, 100000, 2300000)
